Vortices near the Mott phase of a trapped Bose-Einstein condensate 
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We present a theoretical study of vortices within a harmonically trapped Bose-Einstein condensate 
in a rotating optical lattice. We find that proximity to the Mott insulating state dramatically 
effects the vortex structures. To illustrate we give examples in which the vortices: (i) all sit at 
a fixed distance from the center of the trap, forming a ring, or (ii) coalesce at the center of the 
trap, forming a giant vortex. We model time-of-flight expansion to demonstrate the experimental 
observability of our predictions. 
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Atomic clouds in a rotating optical lattice are at the in- 
tellectual intersection of several major paradigms of con- 
densed matter physics. These rotating clouds may dis- 
play a superfluid-insulator quantum phase transition [l[ , 
vortex pinning [2|, frustration ^3j, Josephson junction 
physics and even analogs of the fractional quantum 
Hall effect [f|. Here we explore the theory of vortices in 
such systems, showing how proximity to the Mott insu- 
lator phase impacts the vortex configurations. 

Considering a uniform gas of atoms of mass m in an 
optical lattice rotating with frequency f2, there are three 
macroscopic length scales in the problem: the lattice 
spacing d, the magnetic length I — y/K/mfl, and the 
particle spacing n -1 / 3 , where h = h/2n is Planck's con- 
stant. Even without interactions, the commensurability 
of these lengths leads to nontrivial physics - the single 
particle spectrum, the Hofstadter butterfly, is fractal [6(. 
For interacting bosons, this fractal spectrum leads to a 
modulation of the boundary between superfluid and Mott 
insulating phases @, Further, the vortices in a super- 
fluid on a rotating lattice develop extra structure: their 
cores may fill with the Mott state changing which 
vortex arrangements minimize the energy 0]. 

We consider a harmonically trapped superfluid gas on 
a rotating optical lattice in the single-band tight binding- 
limit close to the Mott state. We choose to study a two- 
dimensional cloud, as it provides the simplest setting for 
investigating vortex physics, and is an experimentally rel- 
evant geometry [10]. The proximity to the Mott state re- 
sults in a nontrivial spatial dependance of the superfluid 
order parameter, and drives a rearrangement of vortices. 

A similar geometry was realized in a recent exper- 



iment [ll|, with the caveat that their shallow optical 
lattice had such a large lattice spacing that they were 
not able to reach the tight binding limit. In principle 
their technique can be refined to explore the physics that 
we describe here. The tight binding limit may also be 
reached through quantum optics techniques which intro- 
duce phases on the hop ping matrix elements for atoms 
in a non-rotating lattice [121 ]. 

In the rotating frame our system is described by the 



rotating Bose-Hubbard hamiltonian [§, [l3| : 
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trix element from site j to site 
tor potential, which gives rise to the Coriolis effect, is 
A(r) — (m/h) (q x fj = ixv (xy — yx), where v is the 
number of circulation quanta per optical-lattice site. The 
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eludes the centripetal potential. In these expressions, /j,q 
is the central chemical potential, u> is the trapping fre- 
quency, fi is the position of site i, d\ (a,) is a bosonic 
creation (annihilation) operator, hi = a}dj is the particle 
number operator for site i, and U is the particle-particle 
interaction strength. The connection between these pa- 
rameters and experiment are given by Jaksch et al. [l3| • 
Here, and in the rest of the paper, we use units where 
the lattice spacing is unity. 

Both the superfluid and Mott insulator can be approx- 
imated by a spatially inhomogeneous Gutzwiller product 
ansatz 13], \^gw) = IlfLi (S„/nl n )0' where i is the 
site index, M is the total number of sites, \n)i is the n- 
particle occupation- number state at site i, and P n is the 
corresponding complex amplitude, with = 1- 

Despite the limitations of being a mean-field theory, the 
Gutzwiller approach compares well with exact methods, 
and strong coupling expansions [l4j . It has also been used 
extensively to understand experimental results [l|, [l5| . 
and is well suited for studying the vortex physics that we 
consider here. 

Using I^gw) as a variational ansatz, we minimize the 
energy with respect to the {/*}. We then extract the 
density p t = ^2 n n|/^| 2 and the condensate order param- 
eter at = (at) = J^n Vn(fn-l)* fn at each site - The 
condensate density — | (a^) | 2 is equal to the superfluid 
density in this model, and is generally not equal to the 
density. 

We use an iterative algorithm to determine the {f„} 
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which minimize the energy. We use a square region with 
L sites per side and hard boundary conditions. We find 
that we must take L much larger than the effective cloud 
diameter so that our solutions do not depend on those 
boundaries. Typically we use 40 < L < 90. For the 
simulation described in figure 1 we impose four-fold rota- 
tional symmetry, but from unconstrained simulations on 
smaller clouds we find that this constraint does not signif- 
icantly change the phenomena. Similar calculations were 
performed by Scarola and Das Sarma [lil ] to analyze the 
case where the single-particle Mott state is surrounded 
by a rotating superfluid ring. 

Since this mean-field theory is highly nonlinear we find 
that the iterative algorithm often converges to different 
solutions depending on the initial state we use. For the 
results shown here we first iterate to self-consistency in 
a parameter region where the solution is unique, then 
slowly change parameter values, using the result from 
the previous parameters as a seed. One should see anal- 
ogous results in an experiment where one adiabatically 
changes parameters. As in such experiments [ljj we ob- 
serve hysteresis. 

We have performed a thorough investigation of a wide 
range of parameters and, as one would expect, we find 
that a basic understanding of the trapped gas can be 
extracted from the phase diagram of the homogeneous 
system, where fl = u) = [Figs. 1(a) and 2(a)]. For 
a sufficiently gentle trap, the gas looks locally homoge- 
neous, and its density at any point r can be approxi- 
mated by that of a uniform system with chemical po- 
tential p(r) = fig — V(r). As a general rule, nontrivial 
vortex structures appear when the LDA superfluid den- 
sity deviates significantly from a typical Thomas-Fermi 
profile. The vortices tend to move to regions where there 
is a local suppression of the superfluid density. 

We illustrate this principle with two examples: in 
Fig. [1] we study the case where the superfluid density 
has a ring-shaped plateau, and in Fig. [2] we consider the 
case where a Mott region sits in the center of the cloud. 

We begin with the nonrotating configuration illus- 
trated by the left half of the figures in Fig. 1 (with 
t/U = 0.06, Ho/U — 0.7). There is a plateau in the super- 
fluid density but not in the total density, and the phase of 
the superfluid order parameter is uniform. Starting from 
this non-rotating configuration we gradually increase the 
rotation speed to v = 0.04, iterating to self-consistency 
at each step. Rather than forming a lattice, the resulting 
vortices form a ring around the central p c peak in Fig. 
1(d). This configuration is favored because it minimizes 
the sum of competing energy costs: the rotation favors 
a uniform distribution of vortices, but the single vortex 
energy is smallest where p c is low. 

As seen in Fig. 1(c), the phase of the superfluid inside 
the ring is essentially constant. This can be understood 
by an analogy with magnetostatics. The velocity v obeys 
an analog of Ampere's law §v-dl = (h/m)N v , where N v , 




FIG. 1: Ring vortex configuration (color online, one- 
column). Comparison between non-rotating (y = 0) and 
rotating (y = 0.04) states of a system characterized by 
(t/U = 0.06, no/U = 0.7). (a) Mean-field phase plot of the 
uniform Bose-Hubbard model. Contours of fixed p and p c , 
are indicated by red and black curves. The superfluid density 
vanishes in the single-particle Mott region labeled "n = 1", 
and increases with lightening shades of purple. The vertical 
orange line represents the LDA parameter-space trajectory 
for the current system, (b) [(d)] Comparison of density [con- 
densate density], (c) Comparison of order parameter com- 
plex phase field. The complex phase is represented by "Hue" . 
Solid and dotted white lines are a guide to the eye. Black cir- 
cles enclose singly-quantized vortices. As seen in (c) and (d), 
vortices form in a circular pattern on the condensate density 
plateau; the density (b) changes only slightly due to rotation. 



the number of vortices enclosed in the contour of integra- 
tion, plays the role of the enclosed current. Neglecting 
the discreteness of the vortices in the ring, the fluid inside 
is motionless, while the fluid outside moves as if all the 
vortices were at the geometric center of the cloud. Even 
with only eight vortices our system appears to approach 
this limit. If one increases the rotation speed, one can 
find a state with several concentric rings of vortices in 
the plateau. Similarly, increasing pa/U can can lead to 
multiple superfluid plateaus, each of which may contain a 
ring of vortices. This structure of nested rings of vortices 
is reminiscent of Onsager and London's original proposal 
of vortex sheets in liquid helium [l8| . 

Our second example of nontrivial vortex structures is 
illustrated in figure [21 where the LDA predicts a super- 
fluid shell surrounding a Mott core. Rather than forming 
a lattice of discrete vortices, one expects that this system 
will form a "giant" vortex [19| when rotated: the vortices 
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FIG. 2: Giant vortex (color online, one-column). Compar- 
ison of a vortex lattice far from the Mott regime (t/U = 0.2) 
and a giant vortex system where the Mott phase occupies 
the center of the cloud (t/U = 0.03). (a) Mean-field phase 
plot for the uniform Bose-Hubbard model with vertical or- 
ange parameter-space trajectory representing a system with 
(t/U — 0.03, (jq/U = 0.3). (b) Comparison of order parame- 
ter complex phase fields, (c) [(d)] Comparison of density and 
condensate density where t/U = 0.03 [t/U = 0.2]. 



occupying the Mott region, leaving a persistent current 
in the supcrfiuid shell. The energy barriers for changing 
vorticity are particularly high, so we generate the rotat- 
ing state in two stages. We start with a non-rotating sys- 
tem (v = 0) at weak coupling (t/U = 0.2, [Iq/U = 0.3), 
gradually increasing the rotation to v — 0.032, where we 
find the square vortex lattice illustrated in Figs. 2(b - 
left) and (d). We then adiabatically reduce t/U from 0.2 
to 0.03. As we reduce t/U, the central p° drops, while p 
approaches unity there. Eventually we see a Mott regime 
at the center of the cloud. During the evolution, we find 
that 8 of the vortices escape from the edge of the trap, 
while four of the vortices coalesce at the center of the trap 
and effectively form a vortex of charge 4. Such a dense 
packing of vorticity would be unstable in the absence of 
the optical lattice. For larger systems with higher rota- 
tion rates one finds giant vortices with larger circulation. 

Due to its multiply connected topology, a ring, such 
as the one formed here, is one of the archetypical geome- 
tries used in theoretical discussions of superfluidity [2(| ■ 
There are several experimental schemes for creating a 
ring-shaped trap [21(, and many theoretical studies of 
giant vortex formation stabilized by a quadratic-plus- 
quartic potential [22] |. Here the multiply connected ge- 



ometry is spontaneously formed by the appearance of the 
Mott state in the center of the cloud. As was found by 
Scarola and Das Sarma [l6j], this Mott region effectively 
pins the vortices to the center. 

By changing t/U one may study a few other interest- 
ing structures. For example, one can engineer a situa- 
tion where a central superfluid region is surrounded by 
a Mott ring followed by a superfluid ring. At appropri- 
ate rotation speeds one produces a configuration which 
has properties of both the states seen in Fig. [T] and in 
Fig. [2 One will find no vortex cores (all of the vorticity 
is confined to the Mott ring), the central region will be 
stationary, and the outer region rotates. 

Another interesting limit is found when one decreases 
the thickness of a Mott/superfiuid region so much that it 
breaks up into a number of discrete islands. Small Mott 
islands act as pinning centers, while small superfluid is- 
lands form an analog of a Josephson junction array [23| . 

Detection. Vortex structures near the Mott limit may 
be hard to detect using in-situ absorption imaging. As 
is exemplified by Fig. 1(b), the vortices do not neces- 
sarily have a great influence on the density of the cloud. 
This is principally because near the Mott boundary the 
superfluid fraction becomes small: even though the su- 
perfluid vanishes in the vortex core, the corresponding 
density may not appreciably change. Two other pieces of 
physics also influence the visibility. First, near the Mott 
boundary one can produce vortices with Mott cores 
Depending on the bulk density, this can lead to vortices 
where there is no density suppression at all, or even a den- 
sity enhancement. Second, the lengthscale of the vortex 
core, the superfluid healing length, varies with U/t. For 
both very large and small U/t the healing length is large, 
while at intermediate couplings it is comparable to the 
lattice spacing, possibly below optical resolution. 

We argue that the vortex structures will be much more 
easily imaged after time-of-flight (TOF) expansion of the 
cloud for time t [3] • The density after TOF expansion is 
made of two pieces - a largely featureless incoherent back- 
ground from the normal component of the gas, and a co- 
herent contribution from the superfluid component. The 
coherent contribution forms a series of Bragg peaks 
where each peak reflects the Fourier transform of the su- 
perfluid order parameter. Neglecting interactions during 
time-of-flight, and using Gaussian initial states at each 
lattice site, we calculate the TOF column density. In the 
long time limit Dj 3> Rtf, where Dj — tit/mX, Rtf is 
the radius of the initial cloud, and A is the size of each ini- 
tial Wannier state, the column density of the expanding 
cloud is 
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FIG. 3: Time-of-flight expansion (color online, one- 
column), (a) Column density p scaled by the central column 
density po as a function of space for the ring configuration vor- 
tex state in Fig. Q] after expanding for time t. Positions are 
measured in terms of scaling parameter Dj = ht/mX, where 
A is the initial extent of the Wannier wavefunction. (b) One 
dimensional cut through center of (a). Note the incoherent 
background between the Bragg peaks, (c) Close-up of the 
central Bragg peak, corresponding to the Fourier transform 
of the superfluid order parameter. The 8 dips in the outer 
crest result from the 8 vortices in the initial state. 

where N and N c are the total number of particles and 
condensed particles, respectively. 

The incoherent contribution (N — N c )p(r,t) is sim- 
ply a Gaussian. This is a consequence of the Gutzwiller 
approximation, which neglects short range correlations. 
Adding these correlations would modify the shape of the 
background, but it will remain smooth. The coherent 
part has much more structure. Figure [3] illustrates the 
density pattern which will be seen if the rotating cloud 
in fig. [T] is allowed to expand. 
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